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Preface 

This  report  will  hopefully  help  in  paving  the  way 
toward  a  useful  algorithm  for  the  estimation  of  launch 
vehicle  parameters.  Although  no  such  algorithm  is  contained 
herein,  I  feel  that  I've  identified  some  major  problems  and 
perhaps  have  recommended  reasonable  approaches  to  solutions 
to  those  problems. 

This  work  could  not  have  become  a  reality  without  the 
typing  efforts  of  Ladonna  Stitzel.  Her  help  is  deeply 
appreciated.  I'd  also  like  to  express  special  thanks  to 
my  advisor,  Dr.  Bill  Wiesel.  A  thesis  effort  which 
doesn't  meet  personal  goals  can  be  frustrating,  but 
his  guidance  and  inspiration  made  this  a  fruitful  learning 
experience. 
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Abstract 


A  seven-state  inverse  covariance  (Bayes)  filter  was 
implemented  to  determine  performance  parameters  of  a 
launch  vehicle.  Data  measurements  were  restricted  to 
azimuth  and  elevation  readings,  typical  of  data  from  an 
infrared  sensor  in  geosynchronous  orbit.  Results  of 
this  study  indicate  that  the  magnitude  of  constant 
acceleration,  assumed  to  act  in  the  direction  of  velocity 
can  be  estimated  using  a  seven-state  filrer  (3  states 
each  for  position  and  velocity,  and  a  seventh  state  for 
acceleration) .  The  system  is  unobservable  for  short  arcs 
of  data  if  only  one  observer  is  available.  The  addition 
of  a  second  observer  can  allow  the  system  to  be  observed. 
An  ad  hoc  fading  -  memory  technique,  in  which  confidence 
in  the  seventh  state  estimate  was  decreased,  proved 
unsuccessful  in  estimating  variable  acceleration  of  a 
launch  vehicle.  Further  attempts  at  estimating  variable 
acceleration  with  an  eight-state  filter  (3  states  each  fo 
position  and  velocity,  and  seventh  and  eighth  states 
involving  engine  exit  velocity  and  propellent  mass  flow 
rate)  were  unsuccessful. 


'  ,T*ic  ~ .“^-vsv '<SW  '7~  ' 


I  Introduction 

Wien  the  increasing  number  of  reconnaissance  and  other 
data-gather j.ng  satellites,  a  growing  source  of  observation 
data  is  becoming  available  which  can  be  used  to  supplement 
or  replace  ground-based  or  aerial  data  for  various  requirements. 
In  some  instances,  this  space  capability  allows  for  the 
gathering  of  data  previously  unavailable  from  the  ground  or 
air.  Measurements  of  elevation  and  azimuth  from  satellite 
infrared  sensors  fall  into  this  category.  The  use  of  these 
data  in  estimating  parameters  of  ballistic  missiles  and 
reentry  vehicles  he s  therefore  been  limited. 

$tui..3s  hav<--  been  conducted  in  the  past  which  concentrate 
on  estimation  of  position,  velocity,  and  performance 
parameters  of  maneuvering  reentry  vehicles.  The  filters 
teveloped  in  these  studies  have  used  ground-based  tracking 
radarrs  as  their  primary  measurement  source.  Typical  data, 
then,  has  included  azimuth,  elevation,  and  range. 

This  paper  presents  the  development  of  an  inverse 
covariance  (Bayes)  filter  updated  by  satellite  data.  The 
satellite  data  consist  of  elevation  and  azimuth  angles  (no 
range)  as  detected  from  an  orbiting  infrared  sensor.  These 
data  would  be  taken  during  the  ascent  of  rockets  or  ballistic 
missiles.  Of  primary  interest  in  this  development  was  the 
ability  of  the  filter  to  estimate  accelerations  of  the  target 
vehicle  in  addition  to  its  position  and  velocity. 
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Assumptions : 

(1)  The  data  satellites  were  assumed  to  be  in  geosynchronous 
equatorial  orbits. 

(2)  Accuracy  of  the  infrared  sensor  was  considered  the 
same  for  determining  both  elevation  and  azimuth  angles. 

(3)  The  first  filter  was  assumed  to  have  seven  states; 
x,  y,  z  components  of  position  and  velocity,  and  an 
acceleration  term  resulting  from  thrust.  This  model 
was  thought  to  be  general  enough  such  that  it  might  also 
be  used  to  estimate  deceleration  of  a  reentry  vehicle  due 
to  drag. 

(4)  The  acceleration  term  was  assumed  to  act  along  the 
velocity  direction  vector  in  order  to  keep  the  model  simple. 

Chapter  II  of  this  paper  presents  the  derivation  of 
the  dynamics  equations  for  the  filter.  This  is  followed 
in  Chapter  III  by  the  development  of  observation  relationships 
between  the  heat  emitting  ballistic  missile  or  reentry  vehicle 
(subsequently  referred  to  as  the  target  vehicle)  and  the 
satellite-based  infrared  sensor.  In  Chapter  IV  the  filter 
equations  are  derived.  In  addition,  Chapter  IV  also  contains 
a  discussion  of  the  problems  encountered  during  implementation 
of  the  seven-state  filter.  A  follow-on  filter  with  eight 
states  is  introduced  and  briefly  analyzed  in  Chapter  V. 

Chapter  VI  contains  conclusions  resulting  from  the  development 
work  and  presents  recommendations  for  further  study. 


2 


II  Derivation  of  Problem  Dynamics 


The  equations  of  motion  for  the  target  vehicle  were 
derived  from  the  general  two-body  equation: 


where  u  =  GM 

AAA 

r  =  xi  +  yj  +  zk  (2-2) 


The  state  model  for  the  filter  initially  had  seven 
states?  three  for  position,  three  for  velocity,  and  one 
for  acceleration  due  to  thrust  or  drag. 

The  state  vector,  x,  therefore  was: 


x  = 


x 

y 

z 

v> 

V 


(2-3) 


The  state  vector  propagates  in  time  according  to  the 
vector  differential  equation: 

x  =  F  (x  ( t)  ,  t)  (2-4) 

In  x,  the  three  position  states  denote  position 
components  in  the  x-f  y-,  and  z-  directions  in. a r coordinate 
frame  with  its  origin  at  the  earth's  center  and  the  z- 
direction  being  north.  The  rates  of  change  of  these  states 
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were  given  by  their  respective  velocities: 


x  = 


V. 


X 


( 2-5a) 


y  =  vy  (2-5b) 

x  =  vz  (2-5 c) 

The  time  rates  of  change  of  the  three  velocity  states 
were  derived  from  the  two-body  equation  and  the  assumed 
representation  for  acceleration  due  to  thrust  (or  similarly, 
deceleration  due  to  drag) .  The  acceleration  will  be  discussed 
further,  but  was  basically  modeled  as  a  vector  acting  in  the 
direction  of  velocity  and  having  constant  magnitude  over 
time  increments  between  observations.  Therefore  acceleration 
was  given  by: 

a  =  a  (2-6) 


Adding  this  to  the  two-body  equation. 


r 


-3  r  +  a 


(2-7) 


yields  equations  for  the  time  rates  of  change  of  the 
velocity  states: 


v 


z 


x  +  a 


y 

—  y  +  a 

r 


z  +  a 


(2-8a) 

(2-8b) 

(2-8c) 
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As  stated  previously,  the  seventh  state,  acceleration, 
was  modelled  as  constant  over  time  intervals  between 
observations.  It  was  assumed  to  act  always  in  the  direction 
of  the  velocity.  The  significance  of  this  implementation 
was  that  this  model  could  be  used  to  represent  acceleration 
during  ascent  (a  positive  value)  or,  with  a  change  of  sign, 
could  represent  deceleration  of  a  reentry  vehicle  descending 
through  the  atmosphere.  This  implementation  would  allow 
the  same  filter  structure  to  be  used  from  lift-off  to 
reentry  vehicle  impact. 

With  acceleration  constant  between  measurement  updates, 
its  time  rate  of  change  was  given  by: 

a  =  0  (2-5) 

The  vector  F  was  formed  as: 


(2-10) 


With  an  assumed  set  of  initial  conditions  (x(t0),tQ), 
the  vector  equation 


X  =  F  (x,t) 


(2-11) 
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can  be  numerically  integrated  to  obtain  a  nominal  trajectory. 
This  trajectory  is  assumed  to  be  known  based  on  assumed 
initial  conditions,  but  the  true  initial  conditions  differ 
slightly  from  those  assumed,  therefore 


x  (t)  =  xQ  (t)  +  <$x  (t) 


(2-12) 


where 


x ( t )  =  the 
Xq (t) =  the 
6x(t)=  the 


assumed  nominal  trajectory 
true  trajectory 

deviation  from  true  trajectory 


Differentiation  gives 


dt  x0(t) 


+  6x(t) 


(2-13) 


which,  together  with  eq  2-11  gives 

ft  "x0(t)  +  ft  6x(t)  <x0(t)  +  6x (t) )  (2-14) 

where  the  time  dependence  of  F  has  been  incorporated  into 
x(t).  Consider  the  first  component  of  eq  2-14: 

dt  xo(t)  +  ft(Sxo^  =  fo{^x0^t^+5  x0  ^  :X1  (t)  +ixi  (r)  xn^ 

+dXn(tj]  (2-15) 

Letting  t  be  fixed  gives  the  right  side  of  eq  2-15  as 

f0(V'0  +  a'  Xl+  6 - Xn  +  t)  (2-16) 
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Application  of  Taylor's  theorem  gives 


3f0  ,  „afo  3f0 


C3*i 

r3xn 

xo 

O 

* 

xo 

xi 

• 

X1 

• 

• 

X1 

• 

• 

• 

X 

• 

X 

• 

X 

n 

n 

n 

+  0 

i2) 

(2-17) 


where  0(2)  represents  terms  of  second  and  higher  orders 
which  will  be  relatively  small  as  long  as  a,  B  /  •  •  •  Y  are  small. 
So,  to  first  order 


dt  x0(t)+dt6x0(t) 


f0 fx0'xl' *  *  *xnJ 


5xo +  *  6xi + 

xo 


(2-lSa) 


*  3Xn! 


In  a  similar  manner,  the  second  component  of  eq  2-15  is 


§t:Vt)+dt6*i(t)  =  fi(*0'*i-“V 


3f  1 

fix.  +  — 


3fl 

fix.  +  ...+  r  rr  • 


(2-18b) 


But 


X(t)  =  F(x  { t )  ) 


(2-11) 


So  recalling  eq  2-14,  eq  2-18  can  be  v/ritten  as 


dt6x0 (fc) 


f  f 
O  Q 
3x0  3xx 


3  1  3f 


..  3f0 


6xq ( t) ' 


6xx(t) 


(2-19) 


St«Vt) 


x0(t) 

X1  it) 


1 6x  (t) 
n 


xn(t) 


So,  to  first  order,  6x(t)  satisfies  the  time-varying 
linear  differential  equation: 


6xU) 


where 


<5x(t)  =  A  [xQ  (t)  6x  {  c) 
A^xQ(t)j  is  defined  by 

L  J  3*.  t  =  T 


(2-20) 


(2-21) 


X  =  X  (t) 


A,  then,  is  a  time-varying  matrix  whose  elements  are  functions 


of  x  and  t.  The  derivation  of  all  elements  of  A  is  contained 


in  Appendix  A 


The  elements  of  the  matrix  A  were  verified  by  conducting 


a  numerical  check  given  by: 


_  F. (x.+6,t)  -  F  <x,t) 

i3  '  -i~;l - J - 4 - 


(2-22) 


where: 


Aij 


th  th 

=  the  element  in  the  i  row  and  j  columns 
=  a  deviation  on  the  order  of  10~^  or  smaller 


F.(x,t)  =  the  itA  component  of  F,  evaluated  at  x 
—  l  th  *— 

Fi(Xj+6/t)=  the  i  component  of  F,  calculated  when  6  has 

th  — * 

been  added  to  the  j  state  of  x 

In  words,  this  check  gives  an  approximation  to  the 
elements  of  A  as  a  result  of  small  changes  in  the  state  vector, 
x.  If  the  matrix.  A,  derived  numerically  agrees  with  A 
evaluated  at  the  state,  x,  this  provides  assurance  that  the 
partial  derivatives  taken  in  deriving  A  are  correct. 

Eq  2-20  is  linear,  therefore  solutions  to  it  can  be 
superimposed.  Also,  a  fundamental  set  of  n  solutions 
(n  vectors,  <f>^)  can  be  constructed  from 


♦i  =A4>.(t0) 


(2-23) 


where  4>i(tQ)  is  a  vector  having  1  as  the  ith  component  and 
zeros  elsewhere 

The  general  solution  to  eq  2-20  is  given  by 


6x(t)  =  6x1(tQ)  4^(t)  +  5x2  (t0H2(t)  +...6xn(t0)<j>n(t)  (2-24) 


which  can  be  assembled  in  matrix  form  as 


6x(t)  -  i(t,tQ) Sx(t0) 


(2-25) 


Since  the  columns,  <$k  of  jHt,tg)  satisfy  eq  2-23,  so 
does  the  matrix  itself.  The  state  transition  matrix,  J>,  is 
given  by 

|(t,t0)  =  A ( t  )i(t,t0)  (2-26) 

and  j&(t,tg)  =  i 

where  I  =  the  identity  matrix 
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The  observation  relationships  for  this  problem  were 
derived  with  both  the  observer  and  target  vehicle  in  the 
same,  x,  y,  z  coordinate  frame.  The  measurements  from  the 
observer  were  the  two  angles,  azimuth  and  elevation  as 
shown  in  Figure  1 . 

Azimuth  was  the  angle  measured  in  radians  (positive 
in  the  clockwise  direction)  between  a  local  vertical,  z‘ 
and  a  local  position  vector,  r'. 

Elevation  was  the  angle,  also  in  radians,  between  the 
negative  of  the  position  vector  of  the  observer,  -R,  and 
the  vector  from  the  observer  to  the  target  vehicle,  p. 

The  azimuth  and  elevation  relationships  were  derived 
from  the  general  equations: 


cos  y  = 


a  ♦  b 
a|  |b| 


(3-1) 


sin  y  -  "X  ~ 

III  lb! 


(3-2) 
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=  position  vector  of  observer 
=  position  vector  of  launch  vehicle 

=  position  vector  of  launch  vehicle  relative  to  observer 

=  position  vector  (in  a  plane  orthogonal  to  R) 
of  launch  vehicle  relative  to  O' 

=  elevation  angle;  the  angle  subtended  by  p  and  -R 

=  azimuth  angle;  the  angle  subtended  by  r'  and 
the  line  segment  from  O'  to  z' 


Figure  1,  Illustration  of  observation  angles 


a  unit  vector  normal  to  R  and  containing  point  P 
a  unit  vector  along  R 


Figure  2.  Geometry  for  elevation  derivation 

Elevation  is  the  angle  between  "p  and  -R,  therefore 
from  Figure  2  and  eq  3-1 

cos  t  ,  ill  =  1Ll-2  ,  I_l_e  ( 

l?l  R  R  R  R 

From  eq  3-3 

|s!  =  |rj  cos  <j>  (; 

Introducing  direction  to  eq  3-4 

s  - - I  r  |  cos  <b  C 

|5| 

From  vector  addition 

s  +  t  =  r  (; 

which,  combined  with  eq  3-5  gives 

~~  ]  r  |  cos  4>  +  t  =  r  (: 
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Therefore 


t  =  r  -  — —  |  r  I  cos  <j> 

iRl 


(3-8) 


Combining  results  of  eqs  3-3  and  3-8  gives 


—  —  R  i  i  R  •  r  —  —  R  •  r 

t  =  r - r  -  =  r  -  R  - 


R  M  r 


R  |  |  r 


(3-9) 


From  Figure  2 


(el)  = 


(3-10) 


p  =  r  -  R 


(3-11) 


Therefore  from  ea  3-9 


sin (el)  = 


r  -  R 

\  1  R  i  2 

Ir  -  r| 


(3-12) 


which  gives 


r  -  R 


el  =  sin  1 


R  •  r 


r  -  R 


(3-13) 


So  the  equation 
az  =  cos  1 


for  the  azimuth  is 


r  -  R 

^R  •  r  \ 

^  ir  i 

A 

•  k 

r  -  R 

{r  •r\ 

ym2J 

(3-16) 


The  sign  of  the  azimuth  angle  is  ambiguous  from  eq  3-16, 
therefore  when  the  equation  was  implemented  in  the  filter,  the 
sign  was  determined  by  looking  at  the  z  component  of  the  cross 
product,  R  x  r.  If  the  z  component  was  negative,  the  azimuth 
was  negative. 

From  the  equations  for  the  observation  values,  elevation 
and  azimuth,  it  is  seen  that  they  are  a  nonlinear  function  of  the 
first  three  elements  of  the  state  vector  as  well  as  the 
position  vector  of  the  observer.  In  reality,  the  observer's 
position  would  be  known.  Thus,  the  two  element  observation 
vector,  ~z  can  be  written  as  a  function  of  the  state  vector: 


z(ti)  =  G  (x(ti),ti)  (3-17) 

If  Xq  is  the  true  state  vector,  a  perfect  observer  would 
produce  exact  measurements : 

^ ti^  =  G(xQ (t±) ,t±)  (3-18) 

writing 

x0(t±)  =  x(t.)  +  6x(t±)  (3-19) 
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“5T 


/*  *-**•« .  , 


^ T;r  -c.  *-r '*  Wj^T5V^f«f»^wti^7-i«( 


the  true  error,  e  (t.),  can  be  written: 

Z  X 


eg(ti)  =  z  ( t  i )  -  z()(ti)=  G  (x  (t±)  ti)-G(x0(ti)  ,tL)  (3-20a) 


=  G (xQ (ti) +6x(ti) ,ti) -G(x0(ti) ,t±) 

(3-20b) 

3G  [ 

(3-20c) 


3x 


dx (ti) 


Since  we  only  have  an  estimate  instead  of  the  true 
state  vector 

x(ti)=  xr(ti)  +6x(ti)  (3-21)  ; 

A 

where  xr(t^)  is  an  estimate  of  the  state  vector  and 
W  is  a  re£erenoe  trajectory. 

=  G(xr(ti)+6x(ti)  yt^^)  “G(xr  (t±)  ,ti) 

(3-22) 

6x(ti)  (3-23) 

Therefore,  the  residual,  r(t^),  is  linearly  related  to  the 
correction  , 5x(t. ),  between  the  reference  state  ,  x  (t. ), 

JL  »  X 

A 

and  the  state  estimate,  x(t^) .  The  linearizations  are 
embodied  in  the  H  matrix  given  by 


The  residual  is  given  by: 


r(ti)  =  z(t±)  -  zr(ti) 


which  is  approximated  by 
3G 


r(t±)  a 


3x 


x  (t.) 
r  i 
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E(t.)  =  ^(V 


—  l 


3x(t.)  _ 


(3-24) 


Since  G(t^)  is  only  a  function  of  the  three  position 
elements  of  the  state  vector,  the  last  four  elements  of 
each  row  of  K  are  zero.  The  derivation  of  the  first  three 
elements  of  both  rows  of  H  is  given  in  Appendix  B. 

In  order  to  make  the  residual  a  function  of  the 
correction  at  tQ  instead  of  t^ ,  the  state  transition  matrix 
is  again  introduced  as  in  eq  2-25  . 

Therefore 


r(t±)  «  H(ti)6ac(ti)  =  H  tQ)  6x  (tQ) 


(3-25) 


The  matrix  product,  H(ti)  $  (t.^,  tQ)  ,  is  designated  T  ( t± )  . 

The  derivatives  in  H  were  verified  in  the  same  manner 
as  those  in  A.  The  verification  resulted  from  looking 
at  the  equation 


=  G  (x^+A)  -  G(x,  t) 


(3-26) 


where  H.  =  the  i  column  of  the  H  matrix 


G(xi+A,t)=  the  evaluation  of  G  with  A  added  to  the  itn  state 

G(x,t)  =  the  evaluation  of  G  with  no  A  added  to  any  of 
the  states 


A  =  a  small  value  of  the  order  10 
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IV  Seven  -  State  Bayes  Filter 

This  chapter  presents  the  derivation  of  equations  for 
the  Bayes  filter  and  the  results  of  their  implementation 
for  the  seven  -  state  filter. 

Derivation  of  Equations  for  Bayes  Filter 

The  problem  dynamics  are  given  by  the  vector  equation 

x  =  F  (x,  t)  (2-11) 

as  derived  in  Chapter  II.  Because  the  dynamics  of  the 
problem  are  well  understood,  small  deviations  of  the  state 
vector  at  any  time  can  be.  expressed  as 

6x(t)  =  i(t,tQ)  6x ( tQ)  (2-25) 

This  expression  is  valid  as  long  as  the  6x's  are  small. 

Eq  2-25  was  used  to  get  an  expression  for  the  residual, 
rft^),  as  a  linear  function  of  the  correction  at  the  epoch 
time,  6x(tQ) 

r(t±)  a  H(ti)£(ti,tQ) ox(t0> 

-  T(t.)ix(t0)  (3-25) 

There  is  an  error,  due  to  imperfect  measurements, 
between  the  true  residual  available  if  the  true  state  vector 
were  known,  and  the  residual  approximation  in  eq  3-25, 
therefore 

r(ti)  =  T(ti) 6x(tQ)+e(ti)  (4-1) 
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Rearranging,  the  error  .is  expressed  as 

e(t±)  =  r(ti)  -T(ti)6x(tQ)  (4-2 

Associated  with  each  observation  vector,  is  a 

covariance  matrix,  containing  information  about  the 
accuracy  of  the  data  instruments.  If  the  measurements 
being  processed  are  mutually  independent,  is  a  diagonal 
matrix.  Using  Gaussian  error  statistics,  the  probability 
of  getting  a  particular  error  vector  can  be  expressed  as 

f(e)=  (2tt)“N/2  |2i^  exp  (-%J)  (4-3 

where  N  =  the  number  of  measurements 

£  =  the  covariance  matrix  associated  with  the  data 

_T 

J  =  e  £  e  and  is  a  scalar 
The  principle  of  maximum  likelihood  is  invoked  by 
maximizing  f.  J  is  a  quadratic  form  and  must  be  minimized 
in  order  that  f  be  a  maximum.  This  is  done  by  solving 

£  =  &  *  0  (4'4 

—  T  T  T 

for  6x(tg).  Using  the  identity  (ab)  =  b  a  and 

dropping  the  time  dependence  for  clarity, 

J  =  (r  -  T6x)T0"1(r  -  T6x) 

— T  -1—  —T  -1  —  — T  T  -1—  — T  T  -1  — 

=  rx£  r  -  r  £  l$x  ~  5x  T  £  r  +  Sx  T  £  T5x  (4- 

g  _<p g  _ip_  — 

Noting  that  —(ax)  =  -^(x  a)=  a  ,  the  partial  derivative 
3x  3x 

3  J 

—  ,  can  be  taken  and  set  equal  to  zero: 

3x 


Vs  .  *  WH  ■»'■ 


|~  =  0-{rT£-1T)T  -  TT£_ 1 r  +  (6xTTT£"1T)T  +  TT£“1T6x 


T  - 1  —  T  - 1—  T  - 1  —  T  - 1  — 

=  -T  2  r  -  T  O  r  +  T  2  I6x  +  T  2  TSx  (4-6) 


-IT  -1 

(Note:  0  ir  symmetric,  therefore  (2  )  =  £  ) 


Setting  —  equal  to  zero  and  dividing  by  two  yields 
3x 


T  - 1—  T  - 1  — 

0  =  -T  g  r  +  f2  T5x 


(4-7) 


which  gives 


—  T  -1  -1  T  - l— 

Sx  =  (T  2  T)  T  £  r 


(4-8) 


T  -1  -1 

where  the  inverse,  (T  Q  T)  ,  exists  as  long  as  the 


states  are  observable  over  the  interval  of  interest. 


To  fine  the  covariance  associated  with  6x,  eq  4-8 


can  be  rewritten  as 


6x  =  W  r  +  e. 


(4-9) 


T  -1  -1  T  -1 

where  W  =  (T  £  T)  T  £ 


Using  this  notation,  the  covariance  of  x  (tf )  , 


Px-(V  =  E  (x  -  xo}  (x  "  VT 


(4-10) 


is  rewritten  as 


s;<V  - E 


r_  —  —  —  t  t 

w|(z  -  zQ)  (z  -  zQ)  W 


(4-11) 


and  since  W  is  deterministic, 


p;<v 


=  W  E 


—  —  —  —  T  T 

(Z  -  Z  Q  )  (Z  “  Z  Q  )  W 


(4-12) 


Since  £  is  the  covariance  of  the  measurements  eq  4-12  can 


be  written  as 


pj(to)  =  £2r 


(4-13) 


■  ■ 


Expanding  and  simplifying; 


T-l  -IT-1  I  T-l  -lT-ll 
l~(tQ)  =  ff  £  XT)  T  Q  XQ  (T  O  T)  T A£  1  I 

T  -1  -IT  T  -1  T  r  T-l  -1  T 

=  (T  0  T)  ‘T  O  )  (T  fi  T)  I 


£x(to)  = 


(4-14) 


These  equations  can  be  altered  somewhat  to  handle  data 
in  a  sequential  manner.  If  new  data,  z,is  to  be  added  to 

A 

an  old  estimate,  x(t0~) ,  and  its  covariance,  P(tQ  ),  the 
old  estimate  can  itself  be  treated  as  data.  The  observation 
relations  for  the  estimate  are 

x  =  JE  x  (ty)  (4-15) 

z  =  G  (x,  t)  (4-16) 

where  I  =  the  identity  matrix 


The  augmented  matrices  are  formed  as 


T 

— aug 


2z 


(4-1?) 


-aug 


£<V  0 


L° 


2Z 


(4-17b) 


'aug 


x 


ref 


(417c) 
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v 

ref  =  the  assumed  reference  state. 

The  '  z'  subscript  denotes  a  vector  or 
matrix  associated  with  the  new  data. 

The  'aug'  subscript  denotes  an  augmented  vector 
or  matrix. 


The  updated  covariance  is  given  by 

+  T  -1  -1 

P  =  (T  0  T  ) 

—  —aug  -^aug  —aug 


-l  T  -l  -l 

=  (P  X  (-)  +  Tlz  0  z  Tz) 


(4-18) 


The  correction  to  the  state  vector  is  found  similarly: 


T  -1  -l  T  -1  — 

(T  0  T  )  T  0  r 
—aug  *  aug  —aug  —aug  —  aug  aug 


(P-1  (-)  +  (P-1(-)  (x~  -  *ref)+ 

(4-19) 


Use  of  Equations  in  Filter  Implementation 

Starting  with  a  good  representation  for  the  state 
vector  as  a  reference  solution,  xr(ti-1),  integrate 

x  =  F(x  (t.  )  t)  (2-11) 

X  X  X  f 

to  the  time  of  the  first  measurement,  t^.  Also  integrate 
£  =  A  £(t,ti<_1)  (2-26) 

to  get  the  state  transition  matrix  from  time  t^_^  to  t.^. 

The  residual,  ,  can  be  obtained  from 

r(ti)  =  z(tA)  -  G(xr (t±) (3-22) 
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Obtain  current  values  for  T(t^)  and  0: 


T(ti)  =  H(ti) 


I(t., 

0  1 


t . 

l- 


(3-25) 


(constant  for  this  implementation) 


■The  covariance  is  updated  by  including  the  information 
from  the  measurement: 


P(t 


+ 

i-1 


) 


(t 


i-1 


) 


Tr(ti)0“1T(ti)| 


-1 


(4-13) 


The  correction  to  the  reference  state  vector  is  given  by 


6x 


X  _ (t .  ,  . 

ref  x-1) 


+ 


(4-19) 


This  correction  is  added  to  the  reference  state,  x  (t.  . 

r  i — J-  /  t 

to  get  a  new  reference  state.  This  procedure  is  continued 
until  the  elements  of  the  state  correction,  6x(t^_1), 
and/or  the  residual,  r(ti),  meet  convergence  criteria. 

When  this  occurs,  the  final  xr(ti„j)  is  integrated  to 
time  t.,  where  it  becomes  the  reference  state  to  be 

i 

integrated  to  the  next  measurement  time.  The  covariance 
is  also  propagated  forward  in  time. 


<V  -  i<Wi>  iT(ti'ti-x) 


(4-20) 


The  iterations  again  proceed,  using  the  next  observation. 
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The  Bayes  filter  algorithm  lends  itself  well  to  the 
performance  parameter  prediction  problem.  It  can  be  started 
with  only  a  guess  for  the  reference  state.  The  inverse  of 
the  covariance  matrix,  ^-1(t~),  can  be  initialized  as  a 
zero  matrix,  indicating  there  is  no  a  priori  knowledge  of 
the  system.  This  is  actually  the  case  when  looking  at 
launch  data  from  a  noncooperative  ballistic  missile. 

With  P  *(t^)  equal  to  zero  the  filter  relies  only  on 
measurement  data  to  make  its  corrections  for  the  first 
set  of  measurements. 

Although  the  Bayes  algorithm  is  sequential,  it  can 
easily  handle  batches  of  data  in  a  least  squares  mode. 

This  flexibility  was  found  to  be  useful  for  handling 
launch  data  and  will  be  discussed  further. 


Initial  Implementation  of  the  Filter 


trajectory 


launch 

vehicle 


observer 


Figure  4.  Initial  positions  of  target  and  observer 


The  observer  in  this  problem  was  assumed  to  be  in  a 
geosynchronous  equatorial  orbit.  Tracking  sequences  were 
initialized  with  the  observer  positioned  on  the  x  axis  as 
shown  in  Figure  4.  This  caused  no  loss  in  generality 
since  the  coordinate  system  chosen  was  arbitrary  and  the 
target  was  initialized  from  a  point  not  also  on  any  of  the 
coordinate  axes.  During  tracking,  the  observer  progressed 
counterclockwise  along  its  orbit. 

Filter  performance  was  checked  by  first  determining 
its  ability  to  determine  zero-acceleration  and  constant 
acceleration  before  looking  at  increasing  acceleration  from 
launch  data.  The  filter  was  able  to  detect  and  correct  a 

-7 

10  perturbation  to  the  x  position  state  on  an  unaccelerated 
trajectory  as  shown  in  Table  1.  The  correction  was  based 
on  8  measurements  at  0.2  time  unit  (approximately  160 
second)  intervals. 

With  the  same  initial  conditions  the  filter  correctly 
detected  a  constant  acceleration  of  10  ^  in  addition  to  a 

-7 

perturbation  to  the  x  position  state  of  10  .  This  correction 

was  based  on  8  measurements  at  0.2  time  unit  (approximately 
160  second)  intervals.  Each  measurement  consisted  of  an 
elevation  and  an  azimuth.  It  should  be  noted  that  for 
each  of  these  cases  the  acceleration  was  constant  which  is 
how  it  was  modeled  in  the  filter.  Also,  the  time  interval 
between  measurements  was  quite  large.  Time  intervals  of 
this  magnitude  are  unacceptable  for  estimating  acceleration 
during  the  ascent  stages  of  a  missile  launch.  The  reason 
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Table  I.  Good  correction  to  first  state 


Filter 

A 

x  initial 

Ax 

from 

true 

Corrections 

1st 

HI 

1111 

HI 

1 

.3026186 

+10'7 

-.100004 

-.697 

eg 

-.603 

-10 

2 

.8637757 

0 

.510483 

-.156 

-10 

.271 

-10 

3 

.3026185 

0 

.152964 

-.475 

-ii 

.218 

-10 

4 

.7316202 

0 

-.732311 

-ii 

-.282 

-10 

-.472 

-10 

5 

.6471292 

0 

'  .168399 

-ii 

.232 

-10 

.449 

-11 

6 

.7316202 

0 

.120171 

-ii 

.819 

-11 

.117 

-10 

7 

.0 

0 

.234957 

-ii 

.120 

-10 

.895 

-11 

Notes:  Corrections  based  on  3  sets  of  8  observations 
Observations  at  0.2  time  unit  (==  160  sec) 
intervals 

Unaccelerated  trajectory 


Table  II.  Good  correction  to  first  and  seventh  states 

Filter 

Ax 

Corrections 

A 

x  initial 

from 

true 

1st 

mm 

3rd 

exp 

1 

.3026186 

+10-7 

-.9999286 

-7 

-.398 

-10 

-.501 

-10 

2 

.8637757 

0 

-.7829917 

-12 

-.926 

-11 

.249 

-10 

3 

.3026185 

0 

-.2544198 

-12 

-.149 

-11 

.  192 

-10 

4 

.7316202 

0 

-.4013456 

-12 

-.550 

-10 

.362 

-10 

5 

.6471292 

0 

-.1409622 

-11 

.297 

-10 

.156 

-11 

6 

.7316202 

0 

-.1307843 

-11 

.135 

-10 

.845 

-11 

7 

.0 

-10“6 

.9999999 

-6 

.210 

-10 

.782 

■ 

Notes:  Corrections  based  on  3  sets  of  8  observations 

Observations  at  0.2  time  unit  (=  160  second) 
intervals  _g 

Trajectory  acceleration  constant  -10 
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for  this  is  that  over  large  time  intervals  on  a  launch, 
the  vehicle  acceleration  would  vary  too  much.  This 
variance  would  make  the  filter's  constant  acceleration 
approximation  less  valid. 

It  was  found  that  with  just  one  observer  the  geometry 
of  the  problem  takes  on  increased  significance.  A 
problem  occurs  when  there  is  too  little  change  in  the 
elevation  and  azimuth  measurements.  This  can  occur  if 
the  time  interval  is  too  small,  if  the  number  of  measure¬ 
ments  is  too  few,  or  if  the  target  vehicle  is  traveling 
away  from  the  observer. 

The  geometry  problem  was  highlighted  during  a  zero- 
acceleration  trajectory.  The  interval  between  measurements 
was  10  seconds  and  measurements  were  processed  in  batches 
of  35.  The  filter  was  unable  to  converge  on  a  solution 
having  zero  acceleration.  Instead,  on  each  subsequent 
correction  the  magnitudes  of  the  corrections  were 
approximately  twice  the  previous  correction  for  each  state 
as  seen  in  Table  III. 

This  suggests  that  the  filter  was  unable  to  distinguish 
the  target's  position  along  the  line  of  sight  vector  as 
shown  in  Figure  5.  The  true  target  position,  r,  was 
indistinguishable  from  alternate  targer  positions,  r|, 
which  would  result  in  the  same  values  for  elevation  and 
azimuth. 

To  deal  with  the  range  observability  problem  a  second 
observer  was  added  to  the  filter  to  give  it  essentially 
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Table  III.  Observability  problem  with  one  observer 


Filter  Ax 

from 

initial  true 


Corrections 


.5168881 

2  -.5144469 

3  .8695028 

4  .4330464 

5  -.4125317 

6  .7156707 

0.0 


0 

0 

0 

0 

0 

0 

+  .001 


-.572  -5 
-.114  -5 
.231  -5 
-.137  -3 
-.168  -4 
.225  -4 
-.947  -3 


Notes:  Corrections  based  on  35  observations 
Observations  at  10  second  intervals 
Zero  acceleration  trajectory 


Table  IV.  Increased  observability 
with  two  observers. 


Filter 

A 

initial 


.5168881 

-.5144469 

.8695028 

.4330464 

-.4123517 

.7156707 

0.0 


Ax 

from 

true 


0 

0 

0 

0 

0 

0 

+  .001 


1st  correction 

exp 

-.1695 

-6 

-.3243 

-6 

.5529 

-7 

.2563 

-5 

.4691 

-5 

-.6703 

-6 

-.9954 

-3 

Notes:  Corrections  based  on  35  observations 
Observations  at  10  second  intervals 
Zero  acceleration  trajectory 


observer  position  vector 
true  target  vehicle  position  vector 
alternate  target  vehicle  position  vectors 
vector  from  observer  to  target  vehicle 


a  "stereo"  view  of  the  target  vehicle's  location.  The 

second  observer  was  assumed  to  be  in  a  geosynchronous 

equitorial  orbit,  but  was  positioned  90  degrees  ahead  of 

the  first  observer  as  shown  in  Figure  6.  . 

With  the  same  initial  conditions  the  trajectory 

was  again  run  with  the  addition  of  the  second  observer. 

The  filter  essentially  corrected  to  the  correct  solution 

in  one  iteration  as  shown  by  Table  IV. 

One  more  trajectory  with  constant  acceleration  was 

tested  before  proceeding  to  launch  trajectories.  This 

final  constant  acceleration  trajectory  started  from 

initial  conditions  similar  to  those  for  a  launch  trajectory. 

Data  was  input  every  second  and  ten  measurements  were 

processed  simultaneously.  The  filter  was  allowed  to  iterate 

using  the  first  batch  of  data  to  see  if  it  would  converge  on 

a  constant  acceleration  solution.  The  first  correction 

was  best,  after  which  subsequent  corrections  shifted  away 

from  the  proper  solution  as  seen  in  Table  V.  Of  particular 

th 

interest  was  the  instability  demonstrated  in  the  10  and 

11  corrections.  These  two  corrections  were  normalized 

and  their  dot  product  was  -0.999991289  which  indicates 

the  vectors  were  parallel.  Increasingly  greater  corrections 

with  opposite  signs  were  being  made  along  the  same  vector 

in  state  space,  characteristic  of  an  oscillating  divergence. 
T  —1 

The  matrix,  T  2  T,  became  noninvertable  following  the 
11  correction.  The  reason  for  this  instability  was 
not  determined. 
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Table  V.  Demonstrated  Filter  Instability  (7-state) 
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Based  on  the  filter's  consistent  ability  to  calculate 
a  good  first  correction  with  constant  acceleration  data, 
filter  performance  with  variable  acceleration  data  was 
investigated.  Data  consisted  of  elevation  and  azimuth 
measurements  as  described  in  Appendix  C  generated  for 
various  time  intervals.  Because  acceleration  was  modeled 
as  constant  in  the  filter,  compensation  had  to  be  made  to 
enable  the  seventh  state  to  change  with  varying  acceleration 
An  ad  hoc  fading  memory  was  added  to  the  filter  in 
which  elements  of  the  covariance  matrix  were  deweighted 
to  reflect  decreased  confidence  in  the  seventh  state. 

This  was  achieved  by  pre-and  post-multiplying  the  inverse 
covariance  matrix,  P  ^(t^),  by  a  diagonal  matrix  just 
after  propagation.  The  elements  of  the  diagonal 
deweighting  matrix  consisted  of  0.99  for  the  first  three 
elements,  0.98  for  the  next  three  elements,  and  the 
seventh  element  taking  on  various  values  less  than  1  to  vary 
filter  performance.  (ref: 8) 

For  filter  attempts  at  estimating  variable 
acceleration,  measurements  were  provided  at  one  second 
intervals  and  were  introduced  to  the  filter  in  sets  of 
five  measurements.  Sets  of  five  were  used  because  of  the 
ease  of  implementation.  As  a  lower  bound,  the  filter  had 
to  have  an  initial  set  of  at  least  four  measurements  in 
order  to  make  all  seven  states  observable.  As  an  upper 
bound,  the  larger  the  number  of  measurements  per  set, 
the  more  the  acceleration  would  vary  over  the  longer  time 


interval.  This  would  make  the  filter's  constant  acceleration 
approximation  less  valid. 

Regardless  of  the  value  selected  for  the  seventh 
element  of  the  deweighting  matrix,  the  filter  estimates 
were  similar.  The  first  five  acceleration  estimates 
increased  linearly  with  values  less  than  the  actual  values 
as  shown  in  Figure  7.  The  filter  would  then  respond  to 
the  increasing  residuals,  overcompensating  with  excessive 
estimates  for  acceleration. 

In  view  of  these  results  an  attempt  was  made  to 
implement  a  better  model  of  acceleration  in  the  filter. 

An  additional  state  was  required  to  do  this.  The  development 
and  results  are  contained  in  the  following  chapter. 
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V  The  Eight  -  State  Filter  Model 


Because  of  the  instability  demonstrated  by  the  seven  - 
state  filter,  the  filter  representation  for  acceleration 
was  changed  to  get  a  better  depiction  of  reality.  Treating 
acceleration  as  constant  between  updates  was  attractive 
because  it  would  have  been  conceptually  general  enough 
to  allow  the  filter  to  process  data  from  the  ascent  of  a 
launch  vehicle  as  well  as  from  the  descent  of  a  reentry 
vehicle.  The  performance  of  such  a  filter  negated  the 
possibility  of  using  sc  general  a  model. 


Derivation  of  Equations 

If  thrust  is  assumed  constant  for  each  stage  of  a 
launch  vehicle  in  accordance  with  usual  theory  (Ref  4:369), 
the  acceleration  of  the  vehicle  can  be  described  by  the 
equation 


a 


T 


V  • 
e 


m 


(mQ-mt) 


(5-1) 


where  aT 
m 


m 


0 


t 


=  acceleration  due  to  the  thrust 
=  the  propellant  mass  flow  rate  (known  as  B  in 

some  literature) 

-  the  original  mass  of  the  missile  (including 
propellant) 

=  time  in  any  consistent  units 
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If  implemented  in  this  manner,  the  filter  would  gain 
three  new  states  -  V  ,  m,  and  ii\q  .  However,  if  the 
numerator  and  denominator  are  divided  by  mg,  aT  becomes: 


am  =  V 
T  e 


(>S j-  t) 

mo 


e  (1-Mt) 


(5-2 


where  M  =  the  relative  mass  flow  rate. 


This  reduces  the  number  of  additional  states  to  two, 
making  an  eight-state  filter  with  the  state  vector: 


The  resulting  equations  of  motion  are: 


x  =  F  (x(t)  ,t)  = 


-x/R  +  aT  vx/v 
-y/R3  +  aT  vy/v 

-z/R-3  +  a_  v  / v 


jpw»w  rv-^-wi  ^m. 


where,  as  before: 


R  =  x2  +  y^  +  z2 

v  =  v2  +  v2  +  v2 
x  y  2 

TT  M 

a  =  V 

T  e  1-Mt 


Further,  the  new  states  result  in  changes  to  the  A 
matrix.  Aside  from  A  now  being  an  8x8  matrix,  there  are 
six  new  non-zero  terms: 
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(5-6) 


v  V, 


(5-7) 


(5-8) 


(5-9) 


A 


6,8 
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M 


(5-10) 


The  added  elements  to  both  G  and  H  due  to  the 
additional  state  are  zeros.  Other  than  modifying  matrix 
dimensions  and  indices  for  compatibility  with  the  eighth 
state,  no  further  changes  to  the  filter  were  required. 


fiiSi 
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Filter  Checkout  and  Performance 

The  filter  was  initially  given  data  with  one  second 

between  observations  (at  each  observation  time  the  filter 

was  given  elevation  and  azimuth  data  from  both  observers) . 

When  these  data  were  processed  to  produce  a  least  squares 

estimate,  both  exit  velocity,  V  ,  and  relative  mass  flow 

rate,  M,were  observable.  However,  when  the  time  interval 

between  observations  was  increased  to  ten  seconds,  the 

T  -1 

resulting  summation  of  T  g  T  was  ill-conditioned  and 

couldn't  be  inverted.  The  diagonal  elements  corresponding 

13 

to  the  position  states  were  of  the  order  10  ;  those  for 

1 8 

velocity  were  of  the  order  10  ;  while  those  for  Vg  and 

2  12 

M  were  10  and  10  ,  respectively. 

As  a  result  of  these  numerical  difficulties,  the 
seventh  and  eighth  states  were  reinspected  for  possible 
changes.  Magnitudes  of  VQ  were  typically  200  to  300 
while  those  for  M  were  approximately  0.005.  Magnitudes 
for  position  and  velocity  were  from  0  to  2.  By  multiplying 
Ve  by  M  in  the  numerator  of  the  acceleration  expression,  the 
term  is  changed  from 

M 

at  e  ( 1+Mt) 


to 


_  V' 
aT  1+Mt 

where  V'  =  V  x  M 
e 


(5-11) 


This  new  parameter  has  dimensions  typical  of  an  acceleration 
and  magnitudes  in  the  range  of  1  -  2 .  When  this  parameter 
was  adopted  into  the  filter  and  the  same  data  was  run  as 
before,  the  estimate  was  again  stopped  by  an  ill-conditioned 
covariance  matrix.  The  diagonal  elements  of  the  last  two 
states  were  both  of  the  order  107,  while  those  corresponding 
to  position  and  velocity  were  as  high  as  before. 

The  same  filter  configuration  was  then  given  data  from 
the  second  stage  of  the  ascent.  Thirteen  observations  were 
processed  simultaneously.  The  covariance  matrix  was  invertible, 
but  the  filter  correction  was  unsatisfactory.  The  filter 
was  not  able  to  identify  and  correct  perturbations  to  the 
initial  conditions  of  the  states.  The  corrections  indicated 
that  the  filter  was  attempting  to  reduce  the  residuals  by 
correcting  all  states  an  equal  amount,  instead  of  just  the 
ones  which  had  been  perturbed  from  nominal  values. 

Based  on  this  result,  observation  data  was  generated 
which  attempted  to  accentuate  the  quadratic  effects  on 
position  due  to  acceleration  in  contrast  to  the  linear 
effects  due  to  velocity.  To  clarify  this  somewhat,  if 
residuals  were  calculated  for  a  trajectory  in  which  the 
initial  position  states  had  been  perturbed,  there  would 
essentially  be  a  constant  displacement  in  the  calculated 
observations.  This  would  result  in  residual  values  which 
were  proportional  to  the  original  position  perturbation: 


where  r  =  residual 

Ap  =  perturbation  in  position 

In  contrast  to  this,  if  the  initial  velocity  states  were 
perturbed,  the  residuals  could  be  expected  to  grow 
linearly  with  time: 

r  *  Av  At  (5-13) 

where  At  =  the  change  in  time 

Av  =  perturbation  in  velocity 

And  finally,  with  perturbation  of  the  initial  acceleration 
state (s),  the  residuals  would  increase  approximately 
quadratically  with  time: 

r  cc  Aa  (At)2  (5-14) 

where  Aa  =  perturbation  in  acceleration 

Therefore,  for  small  trajectory  arcs  (small  At)  the  filter 
could  resolve  the  residuals  by  changing  position,  velocity, 
acceleration,  or  any  combination  of  the  three.  Alternately, 
for  long  arcs  of  data,  acceleration  effects  should  be 
dominant.  With  this  reasoning,  a  trajectory  was  generated 
which  varied  in  acceleration  from  a  value  of  1.09  to  20. 

With  small  corrections  to  either  the  seventh  or  eighth 
state  the  filter  calculated  corrections  for  all  states  as 
seen  in  Tables  VI  and  VII.  Given  the  correct  values  for 
the  two  acceleration  states  and  with  the  fourth  state 
perturbed  the  filter  still  made  inappropriate  corrections 
as  shown  in  Table  VI li  The  reason  for  these  problems  was 


not  determined. 


Table  VI. 

8-staJ 

;e  filter,  7th 

state  perturbed 

Filter 

Ax 

Corrections 

A 

xinitial 

from 

true 

1st 

exp 

2nd 

exp 

1 

1.1045353 

0 

-.39607951 

-6 

-.14425749 

-6 

2 

-.6377037 

0 

-.10541036 

-7 

-.24694020 

-6 

3 

0 

0 

-.53911275 

-7 

-.26563227 

-9 

4 

.3391559 

0 

.27334082 

-6 

.16089711 

-6 

5 

.5874354 

0 

.40171909 

-6 

.69744367 

-6 

6 

.5691713 

0 

.48656250 

-6 

.63308899 

-6 

7 

1.090001 

+  10~6 

-.15398878 

-5 

-.10089342 

-5 

8 

_ 

.00051893 

0 

.10604722 

-9 

1.27149168 

-9 

Notes:  180  observations  processed  simultaneously 

Same  data  for  both  corrections 

10  second  interval  between  observations 
Acceleration  increased  from  1.09  to  20 

Table  VII,  8-state  filter,  8th  state  perturbed 


Filter 

/V 

v 

initial 

Ax 

from 

true 

Corrections 

1st 

exp 

2nd 

exp 

1 

1.1045353 

0 

.44489853 

-3 

.17403788 

-3 

2 

-.6377037 

0 

.13866326 

-4 

.15267806 

-3 

3 

0 

0 

-.27868686 

-5 

.29973877 

-4 

4 

.3391559 

0 

-.46312445 

-3 

-.36587537 

-3 

5 

.5874354 

0 

.52541563’ 

-4 

-.34485975 

-3 

6 

.5691713 

0 

.26822451 

-4 

-.23568659 

-3 

7 

1.09 

0 

.49368101 

-4 

.67655352 

-3 

8 

.00051993 

+  10"6 

-.10231345 

-5 

-.71595424 

-7 

Notes:  180  observations  processed  simultaneously 

Same  data  for  both  corrections 
10  second  interval  between  observations 
Acceleration  increased  from  1.09  to  20 
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Table  VIII.  8-state  filter,  4th  state  perturbed 


Filter 

xinitial 

Ax 

from 

true 

Corrections 

1st 

exp 

2nd 

1 

1.1045353 

0 

-.12170983 

-5 

1.79428858 

2 

-.6377037 

0 

.81175116 

-7 

-.75370992 

3 

0 

0 

-.27062575 

-7 

-.90027939 

4 

.3391569 

+10 

.14454243 

-5 

.22434660 

5 

.5874354 

0 

.34219555 

-6 

.86850929 

6 

.5691713 

0 

.21885328 

-7 

.37698359 

1.09 

0 

-.21938892 

-5 

-.29797983 

.00051893 

0 

.39798295 

-9 

.64701346 

Notes:  180  observations  processed  simultaneously 

Same  data  for  both  corrections 
10  second  interval  between  observations 
Acceleration  increased  from  1.09  to  20 


VI  conclusions  and Recommendations 


Conclusions 

(1)  The  seven-state  (3  states  for  position,  3  states 
for  velocity,  and  1  for  acceleration)  Bayes  filter  can 
estimate  acceleration  of  a  launch  vehicle  if  the  acceleration 
is  constant  (including  zero  acceleration) . 

(2)  Study  results  indicate  that  variable  acceleration 
cannot  be  estimated  with  the  seven-state  Bayes  filter 
employing  a  fading  memory  to  compensate  for  the  constant 
acceleration  modelled  in  the  seventh  state. 

(3)  Results  indicate  that  variable  acceleration 
cannot  be  estimated  with  an  eight-state  Bayes  filter  with 
acceleration  modelled  using  engine  exit  velocity,  launch 
vehicle  mass,  and  propellant  mass  flow  rate. 

(4)  If  measurements  from  only  one  observer  are 
available,  data  arcs  yielding  small  changes  xn  the  measure¬ 
ment  angles  may  result  in  an  unobservable  system.  The 
short  data  arcs  can  be  the  result  of  too  short  a  time 
interval  of  observation  or  problem  geometry. 

Recommendations 

Some  alternate  methods  need  to  be  studied  which 
specifically  take  advantage  of  the  type  of  measurements 
available  in  this  problem.  The  azimuth  and  elevation 
measurements  should  result  in  a  good  position  time  -  history 
for  a  launch  vehicle.  Two  possible  approaches  to  using 
this  information  follow. 
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One  possibility  might  be  to  alter  the  existing 
structure  of  the  filter  algorithms  used  in  this  study. 
Typically,  at  the  initiation  of  a  launch,  the  launch  position 
is  well  known.  The  velocity  starts  at  zero,  but  its 
direction  is  not  well  known  a  priori.  The  acceleration 
can  be  said  to  be  greater  than  one  but  is  otherwise  unknown. 
The  basic  idea  would  be  to  correct  the  states  according  to 
their  expected  contribution  to  the  residuals.  Therefore, 
the  acceleration  state (s)  would  be  corrected  (ignoring 
the  other  six  states)  on  the  initial  filter  iterations (s) . 

The  velocity  and  position  could  be  corrected  subsequently 
(velocity  corrected  before  position) .  If  the  launch 
location  was  assumed  known,  position  wouldn't  be  corrected 
at  all. 

Another  alternative  would  be  to  use  a  smoothing  process. 
Elevation  and  azimuth  data  would  be  used  to  get  a  time-history 
launch  vehicle  position.  A  smoothing  process  (ref:5)  could 
be  used  to  get  a  smooth  curve  which  would  be  differentiated 
analytically  or  numerically  to  get  velocity  and  then 
acceleration.  If  numerical  differentiation  is  used  to 
get  velocity,  the  velocity  curve  would  also  have  to  be 
smoothed  before  differentiating  to  get  acceleration.  In 
addition  the  acceleration  might  also  need  smoothing  to 
get  a  good  representation  of  acceleration  for  all  times 
of  interest. 

Further  analysis  should  be  conducted  to  investigate 
the  causes  of  some  of  the  problems  encountered  during 
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this  study.  With  regard  to  observability  problems,  range 
measurements  might  be  added  to  see  if  the  problems 
disappear.  If  so,  the  deletion  of  range  measurements  and 
addition  of  another  observer  could  be  tried  to  see  if 
similar  results  are  obtained.  If,  instead,  the  addition  of 
range  measurements  doesn't  affect  observability,  further 
analysis  would  have  to  be  conducted  to  look  for  the 
real  cause  of  problems. 

Special  attention  should  be  paid  to  the  time  intervals 
of  observation.  Long  term  propagation  of  the  state  trans¬ 
ition  matrix  may  be  erroneous.  This  could  then  be  the 
source  of  many  problems. 

With  the  seven  state  filter,  some  alternate  fading 
memory  techniques  might  be  attempted.  Also,  the  use  of 
pseudo-noise  was  not  tried  in  this  effort  because  of  the 
desire  to  keep  the  model  simple.  In  view  of  results 
obtained  thus  far,  pseudo-noise  addition  is  a  possibility 
which  might  be  investigated. 
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Appendix  A 


Derivation  of  the  A  Matrix 

The  elements  of  the  A  matrix  are  found  by  taking  the 
gradient  of  the  dynamics  matrix  (vector) ,  F. 

A  =  V? 


where  the  elements  are  given  by 


Using  the  state  equations  given  in  Chapter  II,  the  nonzero 
elements  of  A  are: 


(A- 5) 

(A-6) 

(A-7J 

(A-8) 

(A-9) 

(A-10) 
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where 
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Z  = 


3x 

3v 


3uzx 

5 
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Z  = 


ay 

3v 


3uzy 

5 

r 
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z  = 


3z 


-_U_  +  3pz__ 

3  $ 

r  r 
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66 


A 


67 


=  9vz 

=  -vzvxa 

3Vx 

V3 

=  9vz 

-y  v  a 
=  Jz  y 

xy 

3 

V 

=  9vz 

2 

=-v2a  +  a 

3vz 

v3  v 

=  ^z 

=  Is. 

3 1 

V 

2  ,  2 
+y  +z 

2 ,  2 
+v  +v 

y  z 

and  other  variables  are  as  defined  in  Chapter  II 
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Appendix  B 


Derivation  of  H  Matrix 


The  elements  of  the  H  matrix  are  found  by  taking  the 
partial  derivatives  of  elevation  and  azimuth  with  respect 
to  all  elements  of  the  state  vector.  H  is  therefore  a 


2  by  7  matrix  given  by: 


H  = 


3  G 
3x 


3el  ^el  3_e_l  3_el  iel  jjel.  3el 

3x  3y  3z  3v  3v  3v  3a 

x  y  z 

3az  3az  3az  3az  3az  3az  3az 

3x  3y  3z  3vx  3v^  3vz  3a 


The  last  four  columns  of  the  matrix  are  zero. 


(3-1) 


Partial  ■derivatives  of  elevation: 


el  =  sin  1 


(r 

•R) 

I  -  R  |-i 

i  —  i  +  r 

1  1 R 1 

\R.  \  - 

r-R 


(B-2) 


d 

dx 


(sin  1u)  = 


du 

dx 


(B-3) 


lr'  =  'xW+l2  <B"4) 

where  x,y,z  are  position  components  for  the  target  vehicle 

=  ,/R2+R2+r2  (B-5 ) 

x  y  z 

where  R  ,  R,  R  are  position  components  of  the  observer ♦ 


1 


(B-ll) 


d  -1  .  ~  1 

&  (cos  u)  = 


du 


From  eqs  (B-10)  and  (B-ll] 


u  = 


±i£_lRI  +  7 


R  R 


-R(r  * R)  +  - 


R  R 


(B-12) 


the  numerator  is: 


-R  fcR  +yR  +zR  ) 
z  x  J  y  z 

2  2  2 
R  +R  +R 
x  y  z 


+  z 


Let:  DOT  =  xRx+yR^+zRz 


RSQD  -  R^+Ry+Rz 


So  the  numerator  can  be  written:  z  -  R  DOT 

RSQD 


du 

dx 


2 

u  = 


R  DOT  \ 
z  -  _z _ 

RSQD  / 


x- 


R  DOT 
X 

RSQD 


IZ- 


R  DOT 
z 

RSQD 


(-%) 


12- 


R  DOT 
z 

RSQD 


IX- 


R  DOT 

r-  - 

RSQD 


R  DOT 
x 


z- 


-3/2 


RSQD 

(2)  •( 


R  DOT 

r-  - 

RSQD 


R  DOT 
Z 

RSQD 
2 


x- 


R  DOT' 
x 

RSQD 


,.fx. 

1  \  RSQD  ] 


(  R  DOT 
+  (2)  *  RSQD 


j  \  RSQD  j 


(2)  •  z- 


R-  DOT' 
x 

RSQD 


(B-13 ) 


/-  R  R  \ 
x  z  \ 

RSQD  ) 


-R  R 
x  z 


(B-14) 
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APPENDIX  C 


Description  of  Preliminary  Truth  Model 
for  Data  Generation 

The  preliminary  truth  model  used  in  testing  the  filters 
was  admittedly  simple.  Variable  effects  on  the  trajectory 
caused  by  the  atmosphere  {density  changes  and  winds)  and 
variations  in  the  thrust  were  neglected.  These  variations 
were  never  added  due  to  the  level  of  performance  attained 
by  the  filters. 

The  launch  trajectory  was  elliptical  and  was  generated 

using  modified  two-body  dynamics.  During  the  initial 

portion  of  the  trajectory,  the  effect  of  thrust  was  added 

to  the  equations  of  motion  in  the  form  of  an  acceleration. 

The  seven-element  truth  model  state  vector  consisted  of 

three  states  for  position  x,  y,  and  z;  three  states  for 

velocity  -  v  ,  v  ,  and  v  ;  and  a  seventh  state  for 
x  y  z 

acceleration  due  to  thrust  -  a. 

To  model  acceleration, it  was  assumed  that  thrust  was 
constant  for  each  stage  and  that  the  initial 
thrust  to  weight  ratio  was  known  (ref  4:369).  The 
exhaust  exit  velocity  can  be  determined  from 


V  = 
e 


(C-l) 


where  Ve 


=  exit  velocity  of  rocket  engine  exhaust 
(f t/sec) 

=  specific  impulse  of  rocket  engine  (sec) 

2 

=  gravity  (ft/sec  ) 
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The  trajectory  resulted  from  integrating  the  modified 
two-body  equations  of  motion: 


F  =  I  X 


(C-6 ) 


-x/R  +  (a)  (vx) /v 
-y/R3  +  (a) (vy)/v 
-z/R3  +  (a) (v  ) /v 
VeM2/ (1-Mt) 2 


where 

n3  .  2  A  2  ^  2,3/2 

R  =  (x  +  y  +  z  )  ' 

v  =  (v2  +  v2  +  v2) 1/2 

Ve  =  exit  velocity  of  propellant  from  rocket  engine 
M  =  relative  mass  flow  rate,  the  mass  flow  rate  of 
propellant  from  the  engine  divided  by  the  total 
initial  mass  of  the  launch  vehicle,  3/mQ 
t  =  time 

a  =  acceleration  due  to  thrust;  assumed  to  act 
in  the  direction  of  the  velocity  vector 
Initial  conditions  on  many  of  the  states  were  critical 
Criteria  for  the  position  states  were: 

(1)  that  the  square  root  of  the  sum  of  the  squares 

of  the  position  components  equal  1. (distance  unit) 
In  other  words,  the  launch  point  was  on  the 
earth's  surface. 

(2)  Coordinates  were  restricted  from  any  of  the  axes 
in  order  to  keep  away  from  special  cases. 
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(3)  The  launch  point  should  be  xn  the  field  of  view  of 
both  observer  satellites. 

The  criteria  on  the  velocities  were  most  critical: 

(1)  The  magnitudes  had  to  be  small  to  simulate 
near-zero  velocity  at  lift-off,  but  couldn't  be 
zero  because  of  the  need  to  provide  a  direction 
for  acceleration. 

(2)  The  relative  magnitudes  between  velocity  components 
needed  to  follow  closely  magnitudes  between 
position  components  for  a  near-vertical  takeoff. 

(3)  The  velocity  had  to  deviate  slightly  from 
vertical  to  result  in  a  gravity  turn  within  the 
field  of  view  of  the  observers. 

Trial  and  error  was  used  to  come  up  with  proper  initial 

velocity  conditions.  Velocity  magnitudes  approximately 
-3 

10  that  of  the  position  vector  proved  successful.  A 
velocity  magnitude  decrease  of  approximately  one  percent 
from  vertical  was  used  to  result  in  a  gravity  turn  in  the 
direction  of  decreased  magnitude. 

Values  for  the  thrust  profile  were  derived  from 
parameters  for  a  Titan  IIIB  rocket  (ref  6:5-141). 

1st  stage  V^  (exit  velocity)  -  8243.2  ft/sec 

Thrust  -  464,900  lb 

Propellant  mass  flow  r^fe  -  56.398  Ib-sec/ft 

2 

Initial  mass  -  11,275  lb-sec  /ft 
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2nd  stage  VQ  -  10,220.3  ft/sec 
Thrust  -  102,300  lb 

Propellant  mass  flow  rate  -  10.0095  lb-sec/ft 

2 

Initial  mass  -  2735  lb-sec  /ft 


3rd  stage  V  -  9402.4  ft/sec 
Thrust  -  16,000  lb 

Propellant  mass  flow  rate  -  1.7017  lb-sec/ft 

2 

Initial  mass  -  455  lb-sec  /ft 


To  simulate  staging,  time  dependent  conditional 
commands  were  used  to  select  a  different  subroutine  for 
each  stage.  At  the  beginning  of  each  stage  the  acceleration 
would  begin  with  a  nominal  value,  V£M.  Acceleration 
increased,  reaching  its  maximum  just  before  staging  as 
shown  in  Figure  C-l. 
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